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HIGH  ORDER  FINITE  DIFFERENCE  METHODS,  MULTIDIMENSIONAL  LINEAR 
PROBLEMS  AND  CURVILINEAR  COORDINATES 

JAN  NORDSTROM*  AND  MARK  H.  CARPENTERt 


Abstract.  Boundary  and  interface  conditions  are  derived  for  high  order  finite  difference  methods  applied 
to  multidimensional  linear  problems  in  curvilinear  coordinates.  The  boundary  and  interface  conditions  lead 
to  conservative  schemes  and  strict  and  strong  stability  provided  that  certain  metric  conditions  are  met. 

Key  words,  high-order  finite-difference,  numerical  stability,  interface  conditions,  summation-by-parts, 
variable  coefficient 


Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction.  Phenomena  that  require  an  accurate  description  of  high  frequency  variation  in  space 
for  long  times  occur  in  many  important  applications  such  as  electromagnetics,  acoustics  (all  cases  of  wave 
propagation),  and  direct  simulation  of  turbulent  and  transitional  flow;  see  for  example  [l]-[6].  Strictly  stable 
high  order  finite  difference  methods  are  well  suited  for  these  types  of  problems  (see  [7]- [16])  because  they 
guarantee  accurate  results  with  bounded  error  growth  in  time  for  realistic  meshes. 

Most  of  the  development  for  these  types  of  methods  has  considered  constant  coefficient  problems  on  a 
Cartesian  mesh.  In  [17].  [18]  stable  and  conservative  boundary  and  interface  conditions,  were  derived  for 
the  one-dimensional  constant  coefficient  Euler  and  Navier-Stokes  equations  on  multiple  domains.  A  similar 
technique  was  used  in  [19].  [20].  and  [21]  for  Chebyshev  spectral  methods. 

In  this  paper  we  extend  the  constant  coefficient  analysis  in  [17].  [18]  to  scalar  multidimensional  linear 
problems  in  curvilinear  coordinates  including  block  interfaces.  Related  previous  work  includes  investigations 
of  the  metric  derivatives  in  non-smooth  meshes  (see  [22].  [23])  and  the  treatment  of  parabolic  and  hyperbolic 
systems  in  curvilinear  coordinates  on  a  single  domain  [14]. 

The  rest  of  this  paper  will  proceed  as  follows.  Section  2,  will  give  some  basic  definitions.  Section  3 
presents  the  ID  difference  operators  that  form  the  basis  of  the  multidimensional  approximation  treatment. 
Section  4  defines  the  linear  model  problem  and  discusses  well-posedness.  Section  5  provides  an  investigation 
of  the  discrete  problem.  Section  6  illustrates  numerical  experiments  and  in  Section  7  we  summarize  and 
draw  conclusions. 


2.  Definitions.  Consider  the  linear  initial  boundary  value  problem 


wt  =  P(x.  t)w  +  SF(x.t) 

,x  eQ 

w  ~  Sf(x) 

.  x  £ 

Lcw  =  8g(t) 

.  x  g  r 

A  >  0, 
A  =  0. 
A>  0. 


where  P  is  the  differential  operator.  Lc  is  the  boundary  operator.  Q  is  the  domain  and  T  is  the  domain 
boundary.  The  forcing  function  SF:  the  initial  function  Sf.  and  the  boundary  data  Sg  are  the  data  of  the 
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problem;  w  denotes  the  difference  between  a  solution  with  data  /.  F .  g  and  one  with  data  /  -f  Sf ,  F  +  £F. 

g  +  Sg* 

The  semi-discrete  version  of  (2.1)  is 


(wj  )t  “  Q(xj  t  t)wj  d-  (0 

.  xj  e 

A  >  0; 

J 

II 

C^> 

A  =  0. 

LDWj  -  Sg{t) 

GT 

A  >  0. 

where  Q  is  the  difference  operator  approximating  the  differential  operator  P .  SFj  is  the  forcing  function.  Sfj 
the  initial  function.  Ljy  the  discrete  boundary  operator  where  numerical  boundary  conditions  are  included, 
and  Sg  the  boundary  data.  It  is  assumed  that  (2.2)  is  a  consistent  approximation  of  (2.1). 

2.1.  Well-posedness  and  stability.  There  are  many  concepts  of  well-posedness  and  stability;  see 
[24].  Here  we  consider  the  following  definitions. 

Definition  1.  Problem  (2.1)  is  strongly  well  posed  if  the  solution  w  is  unique,  exists,  and  satisfies 

(2.3)  IM&  +  r  IH£*  <  +  Ailin’  +  ||<m£)a}; 

Jo  Jo 

where  Kc  and  r}c  may  not  depend  on  SF,  Sf,  Sg,  and  ||  *  ||n  and  ||  •  ||r  are  suitable  continuous  norms. 
Definition  2.  Problem  (2.2)  is  strongly  stable  if  for  a  sufficiently  fine  mesh  the  solution  wj  satisfies 

(2.4)  \\w\\l  +  f  |Hft A  <  Kje^mWl  +  f  (\\SF\\2n  +  \\Sg\\frdt}.. 

Jo  Jo 

where  K&  and  rjd  may  not  depend  on  SFj,  Sfj,  Sg,  and  ||  •  || a  and  ||  •  ||r  are  suitable  discrete  norms. 

Definition  3.  The  approximation  (2.2)  of  (2.1)  is  strictly  stable  if  the  analytical  and  discrete  growth 
rates  ( see  (2.3)  and  (2.4))  satisfy 

(2.5)  7]d<7]c  +  0{Ax): 


where  Ax  is  the  mesh  size . 

2.2.  Linear  algebra  relations.  For  later  reference  we  define  some  useful  matrix  operations;  see  [25]. 
Definition  4.  Let  A  be  a  p  x  q  matrix,  and  B  be  an  m  x  n  matrix;  then 


A(&  B  — 


(  ao.oF  •••  ao.g-i^  ^ 

\  Gp-i.oB  •  •  •  ap-i;g-iB  J 


The  p  x  q  block  matrix  A  &  B  is  called  a  Kronecker  product.  There  are  a  number  of  rules  for  Kronecker 
products;  see  [25].  In  this  paper  we  will  make  frequent  use  of 


(2.6) 


[A  ®  B){C  <2  D)  =  {AC)  ®{BD), 
{A®B)T  =at®bt, 

{A  <8  B)~l  =  A"1  0  B”1. 


Consider  the  following  matrices. 

(2.7)  B  =  BT  >  0;  C  =  CT  >  0.  D  =  diag(di)  >  0: 
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where  B.  C .  and  B  <0  C  have  the  structure 

’  Bl  1  [  CL 

1  0 

(2.8)  B=  ,C  = 

0  ,  0 

BR  _ 

We  will  need  the  following  lemmas. 

Lemma  1.  Let  C  and  D  be  the  M  x  M  matrices  defined  in  (2.7)-(2.8).  If  the  blocks  Cl-.Cr  have 
dimension  r  x  r  and  the  first  and  last  r  components  in  D  are  constant,  then 

(2.9)  CD  =  DC  =  C1/2DC1/2  >  0. 

Proof:  A  direct  matrix  multiplication  leads  to  (2.9).  □ 

Lemma  2.  Let  the  first  and  last  r  components  Dk  =  diag(di).  k  =  1. ....  N  be  constants;  let  B  have  q  x  q 
blocks  Bl  .  Br  and  C  have  rxr  blocks  Cl-  Cr;  see  (2. 7)  -(2.8).  With  A  =  B&C  and  D  =  diag(Dk)  we  have 


(2.10)  AD  =  DA  =  A1/2DA1/2  >  0 

if  the  first  (Di . ....  Dq)  and  last  (DN_^q_1y  ...,DN)  q-blocks  in  D  are  equal.  Proof:  By  introducing  DL  = 


(Iq  <g  Dx)  with  I  the  identity  matrix,  the  relations  (2.6)  and  lemma  1.  the  upper  left  corner  of  AD  becomes 
(Bl  0  C)Dl  =Bl  0  CDX  =  BL®  DXC  =  DL(BL  ®  C) 

and 

(Bl  0  C)Dl  =  ( B lL/2  g  C1/2)(Bl/2  ®  DXC1/2)  =  (B1/2  0  C1/2)DL(B1L/2  0  C1/2). 

Positive  definiteness  follows  directly  from  the  fact  that  Di  >0.  □ 

3.  The  ID  difference  operators.  The  ID  difference  operators  that  form  the  basis  for  the  multidi¬ 
mensional  difference  approximations  will  be  presented  below,  for  more  information  see  [17].  [18]. 

3.1.  The  discrete  differentiation  operator.  Let  U ,  VU  be  the  numerical  approximations  of  the 
scalar  quantities  u  and  ux  respectively.  The  approximation  VU  of  the  first  derivative 

(3.1)  VU  =  P^QU:  ux-  P~1Qu  =  Tei  \Te\  =  0(Axm,Axn) 

where  \Te\  =  0( Ax™.Axn)  means  that  the  approximation  of  the  differential  operator  is  accurate  to  order 
m  in  the  interior  of  the  domain  and  to  order  n  at  the  boundary.  Typically  we  have  n  —  m  -  1.  The 
summation-by-parts  (SBP)  operator  VU  satisfies 

(3.2)  {U,VV)P  =  UNVN  -  UoV0  -  (: VU,V)P 

where 

(3.3)  (U.V)p  =  UT PV:  P  =  PT,  Q  +  Qt  =  D,  D  =  diag[- 1,0,  ....0,1] 

and  0  <  pmi„AxI  <  P  <  pmaxAxI.  Operators  of  the  SBP  type  arise  naturally  with  centered  difference 
approximations;  for  examples  see  [9];  [15].  [12].  [26]. 
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In  this  paper  we  will  apply  the  first  derivative  operator  twice  to  obtain  the  second  derivative:  i.e..  we 
will  use 

(3.4)  V2U  =  V{VU):  uxx-P~1QP~1Qu=Te:  Te  =  0{Axm.Axp) 

despite  the  fact  that  we  lose  two  orders  of  magnitude  in  accuracy  p  —  m  —  2  at  the  boundaries.  The 
second  derivative  defined  in  (3.4)  satisfies  (3.2)  with  V  =  BVU:  which  is  completely  similar  to  (u;  (bux)x)  = 
ubux\l  -  ( ux .  bux)  obtained  in  the  continuous  case.  For  another  type  of  second  derivative,  see  [17].  [18]. 

3.2.  The  discrete  integration  operator.  The  matrix  P  in  the  SBP  derivative  operator  is  a  discrete 
integration  operator. 

Theorem  1.  Let  the  difference  operator  V  =  P~lQ  defined  in  (3A)-(3.3)  exists  on  the  interval  —1  < 
x  <  1.  Then,  the  matrix  P  is  an  integration  operator  which  satisfies 

(3.5)  j 1  (uv)„dx  =  UTPVV  +  (VUfPV 

where  U.  V  are  the  projections  of  the  continuous  functions  u.  v  onto  the  grid .  Proof:  UT PW+(VU)T PV  = 
(UV)W  =  («w)lti  =  D 

It  is  also  possible  to  prove  the  following  theorem. 

Theorem  2.  Let  the  difference  operator  P~~1Q  defined  in  (3.1)-(3.3)  exist  on  the  interval  —  1  <  x  <  1 
and  be  accurate  to  order  m.  Then,  the  matrix  P  is  an  integration  operator  which  satisfies 

(3.6)  J1  uv  dx=UTPV  +  0(Axm). 

where  U.V  are  the  projections  of  the  continuous  functions  u,  v  onto  the  grid.  Sketch  of  Proof:  The  proof 
has  two  parts.  The  first  part  shows  that  (3.6)  holds  for  general  polynomials.  Next.  Weistrass5  interpolation 
theorem  (see  [28])  is  used  to  show  that  it  also  holds  for  continuous  functions. 

4.  The  continuous  problem.  The  definition  and  specification  of  the  continuous  problem  is  done 
with  Cartesian  coordinates.  After  the  transformation  to  curvilinear  coordinates  we  check  that  the  essential 
mathematical  properties  are  preserved. 

4.1.  Cartesian  coordinates.  The  two-dimensional  (2D)  linear  problem  considered  in  this  paper  is 

Ut  +  Fx  -f  Gy  —  h.  [x,  y]  E  Q;  t  >  0 
(4.1)  u  =  /.  [x,  y]  (E  t  —  0 

Lu  =  g,  [x,y]eSQ,  t>  0. 

where  h.f.  g  are  the  data  of  the  problem.  L  is  the  boundary  operator,  and 

F  =  FT  +  Fvt  F7  =  giu;  Fv  =  -(6ntix  +  ii2«y), 

G  =  G 1  +  Gv ,  Gj  =  a^U.  Gv  =  —(^21^0:  -f  ^22uy)‘ 

The  coefficients  ai.bij  are  known  functions  of  x.y,  and  t.  For  simplicity  we  have  chosen  Q  =  [x.y]  € 
[—1.1]  x  [0;  1];  see  figure  4.1.  For  future  reference  we  also  introduce 

(4.3)  a  =  (ai.a2).  F  =  (F.G).  n  =  (ni.n2 ): 

where  n  is  the  outward  pointing  unit  normal  on  JQ. 

Equation  (4.1)  can  be  thought  of  as  a  model  for  the  Euler.  N a vier- Stokes,  or  Maxwell's  equations. 


NORTH 


qi  -  1 

WEST 

y  — 1 

EAST 


x  =  —  1 


X  =  -f  1 


SOUTH 


FlG.  4.1.  The  computational  domain . 


4.1.1.  The  energy  estimate.  Well-posedness  of  (4.1)  requires  that  we  can  obtain  an  energy  estimate. 
Let 

(4.4)  (u.*;)—  J  j  uv  dxdy .  (u.  u)  —  ||u||2. 

(4.5)  (u,v)e.w  =  /  (u?;)£My<fy.  (u;ti)jE:iy  =  IMli.w 

Jo 

(4.6)  (u.v)N,s  =  J  (uv)N:Sdx.  ||u||^tS  =  (u.u)j\r,s|Mlw,s  =  {u;u)n.s 

denote  the  Li  scalar  product,  the  Li  norm,  the  boundary  scalar  products,,  and  the  boundary  norms  respec¬ 
tively.  The  subscripts  E.W.N.S  refer  to  the  EAST.  WEST.  NORTH,  and  SOUTH  boundaries  (see  figure 
4.1). 

The  energy  method  applied  to  (4.1)  leads  to 

IMIt2  =  -  [(«';  FJ  +  2 FV)E  -  (u,  +  2 F*_ M 

EAST-WEST 

-  [{u.G1  +  2 Gv)n  -  (ti,  G1  +  2Gv)s] 

- - V - ' 

NORTH-SOUTH 

~  [(u;  Fj)  -  (ug;Ff)  +  (u,  Gy)  -  {uy ,  G7)]  +  J(«.  h)  +  (h,  u)] 

GROWTH  1  "  GROWTH  2 

(4.7)  -[(ug;Fv)  +  (Fv;^)  +  (u!/;Gv)  +  (Gv;%)]/ 

DISSIPATION 

GROWTH  1  (GR1)  and  GROWTH  2  (GR2)  in  (4.7)  will  lead  to  a  growth  or  decay  in  ||u||2.  but  will  not 
affect  well-posedness.  Note  that  for  constant  coefficient  problems  GR1  is  zero.  To  bound  |Iu||2  in  time,  the 
first  two  terms  must  be  bounded  using  the  correct  boundary  conditions  and  the  DISSIPATION  (DI)  must 
have  the  right  sign. 

4.1.2.  The  dissipation.  For  a  correct  sign  of  DI  (see  (4.2)  and  (4.7))  the  eigenvalues  of  B  -\-BT  where 


(4.8) 


B  = 


hi  &12 
621  &22 
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must  be  positive;  therefore. 


b\\  >  0.  &22  >  0.  &11&2 


'bu  +  hi 
2 


is  a  requirement  for  well-posedness. 

4.1.3.  Boundary  conditions.  By  integrating  the  “cross  derivative”  components  of  the  first  two  terms 
in  (4.7)  one  obtains 

-  £  u(FJ  +  2Fv)\±\dy-  J  u(Gt  +  2Gv)\ldx  =  (b12  +  621)«2|111|J 
-J  u{(ai  +  (bi2)y)u  -  2buux)\t\dy  -  J  u((a2  +  (&2i)*)«  -  2b22uy)\\dx . 

The  term  (&12  +  &2i)^2|tilo  involves  the  point  values  in  the  four  corners  of  the  computational  domain.  These 
point  values  cannot  be  estimated  (unless  they  are  artificially  specified).  Boundary  conditions  where  the 
viscous  fluxes  ( FV:GV )  are  specified  avoid  that  difficulty  and  lead  to  an  energy  estimate;  therefore  one 
should  specify  the  total  flux  (F  =  FJ  +  Fv .  G  =  G1  +  Gv)  or  the  viscous  flux  at  the  boundaries. 

Consider  the  first  two  terms  in  (4.7)  and  recall  the  definition  (4.3).  At  x  =  ±1  with  n  =  (dbl.  0)  we  have 
the  boundary  conditions 


— ai(— 1 

=  a 

•  n  <  0 

F 

•  n  — 

F 

=  Fw(y-.t): 

-.y-t) 

—  a 

•  n  >  0 

Fv 

•  n  — 

Fv 

=  F^(y.t). 

+ai(+l 

-.y-.i) 

=  a 

•  n  >  0 

Fv 

•  n  — 

Fv 

= 

+ai(+l 

-.y-.t) 

=  a 

•  n  <  0 

F 

•  n  = 

F 

=  Fb(2/:<); 

-a2{x, 

o;t) 

—  a  • 

n  <  0 

F  ■ 

n  = 

G  = 

Gs(M): 

-a2(x. 

o.t) 

=  a  - 

ft  >  0 

Fv  ■ 

n  — 

Gv  = 

G£(M): 

+a2(x. 

1  ,*) 

—  a  • 

n  >  0 

Fv  ■ 

n  = 

Gv  = 

+a2(x. 

1,*) 

=  a  * 

n  <  0 

F  ■ 

n  = 

G  = 

Gjv(<M): 

are  used  at  y  =  0. 1  where  n  =  (0.q=l)  .  The  boundary  conditions  (4.10).  (4.11)  can  also  be  formulated  in 
the  following  more  general  way. 

a  •  n  <  0  =>  F  •  n  =  Fsn  •  n 


a  -  n  <  0  =>  F  •  n  =  Fsn  *  n 

(412)  a •  n > 0  =>  Fv-H.=  Pjf^n 

Boundary  conditions  of  the  type  (4.10).  (4.11).  (4.12)  have  been  derived  in  [29]  and  [20]  for  the  Navier-Stokes 
equations. 

Let  us  consider  the  EAST  boundary  in  detail,  and  let  us  assume  that  a\  is  positive  at  y  =  0.  becomes 
negative  at  y  =  yo  and  remains  negative  until  y  =  1.  Inserting  the  boundary  conditions  (4.10)  into  (4.7) 
yields 

/»1  ry  0 

-  u(FJ  +  2Fv)\^dy  =  -  /  \ai{l.y.t)\u2  +  2uF%dy 

Jo  Jo 

rV  0 

+  /  ai(l.?/.i)u2  -  2uFEdy  =  -  /  |«i |«2  +  2uF%dy 

Jyo  Jo 

(4.13)  —[  |ai|u2  +  2uFE<Iy  —  —  f  \cti\u2  +  2u(cr1F^  -f  Fg)dy. 

J  yo  J  0 
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where  cri  =  ( 1  —  |ai|/ai)/2.  The  requirement  that  ai  changes  sign  in  the  manner  described  above  can  be 
relaxed,  so  (4.13)  is  a  generally  valid  formula.  An  entirely  similar  procedure  at  the  other  (WEST.  NORTH,. 
SOUTH)  boundaries  yields  the  final  result: 

— (u.  F1  +  2 Fv)e  =  — («■  |ai|w)^  “  2{u,  Fe)e: 

(u.  F1  +  2 Fv)w  —  — («.  WiW)w  -  2{u.  Fw)w: 

—  (u.Gt  +  2 Gv)n  =  —  (u*  |«2 1 u)n  —  2(u.  G;v);v: 

(4.14)  («.  GJ  4-  2GV)S  =  -(«;  |fl2|«)5  -  2(u.  65)5; 

where 

Fe  =  ctiFe  +  (1  -  ^1)^ ;  crx  =  (1  —  |ai|/ai)/2. 

Ftv  —  +  (1  —  03)^VV:  ^3  =  — (1  +  |^i  |/ ai)/2r 

Gjv  =  o^Gn  +  (1  —  ^5  )G^f.  cr5  =  (1  —  |a2|/ 0-2)12. 

(4.15)  Gs  =  o-^Gs  4-  (1  —  ^7)^5  .  (T7  =  —  (1  -f  \a2\/a2)/2. 


Inserting  the  relations  (4.14)  and  (4.15)  into  (4.7)  leads  to 


(4.16) 


||«||2<  Y1  -11^/11/  +  GR1  +  GR2  +  DL 

I=E;W:N}S  rfI 


where 


i)E 


Jolaih2^. 

fo^dy  U=1 


T}W  = 


Jo  \ai\u2dy 
fo  u*dy 


fhW^dx  f^ld^dx 

m=  pS  ^  "s~~i vTu 

The  parameters  r}E:  :r}N>  Tjs  are  strictly  positive  if  a\.  a2  are  zero  for  a  finite  number  of  points. 

Time-integration  of  (4.16)  leads  to  an  energy  estimate  of  the  form  (2.3)  if  (4.9)  holds.  Provided  that 
a  solution  exists  (can  be  shown  by  using  the  Laplace-transform  technique  or  via  difference  approximations; 
see  [30]  and  [31]).  we  can  conclude  that  the  following  theorem  holds. 

Theorem  3.  Problem  (4-1),  (4-10).  (4. 11)  is  strongly  well  posed . 

4.2.  Curvilinear  coordinates.  In  this  Section  we  consider  problem  (4.1)  on  a  curvilinear  domain.  By 
introducing  the  transformation  t  =  r.x  =  x{^.rj.r).y  =  y(£.  r}.  r)  and  it's  inverse 

(4.17)  r-t,  r]  =  r}(x.y.t). 

we  obtain  the  transformed  equation 

(4.18)  ( Jn)T  +  ( J{itu  +  £XF  +  *„G))€  +  ( J(r}tu  +  r)xF  +  r}yG))v  =  Jh  +  RHS: 
where 


RHS  =  u(Jt  +  (J£t)t  4-  4*  4“  4-  G({J£y)z  4-  («/%)?/)  • 
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Using  the  metric  relations  causes  the  term  RHS  to  vanish: 

(4.19) 


Jit 

—  ( X1)Vt  —  XtVt]) 

Jrjt  = 

[y^r  - 

-  xiVr) 

Ji* 

=  Vv 

Jiy  = 

%rj 

J  T)x 

=  -Vt 

II 

e? 

*-> 

xt 

J 

=  xzyv-xnys 

j  = 

(Zxly  ' 

~  iyV%) 

In  this  paper  we  will  consider  the  steady  version  of  (4.17);  i.e..  £  =  i(x.y).r)  =  r](x.y).  The  new 
transformed  problem  becomes 


(4.20) 


JuT  +  (F)(  +  (G)r,  =  h:  r>  0; 

u  =  f,  [£.7?]  eft.  t  =  0. 

Lu  =  g ■  [£,»?]  €<5^:  r  >  0.. 


where  h  =  Jh.  f.g  are  the  data  of  the  problem  and  ft  =  [£.  rj]  £  [-1. 1]  x  [0. 1].  The  new  transformed  fluxes 
are 


(4.21) 
where 

(4.22) 


F=  J{F-V£)  =  FI  +  FV;  FJ-a  lU;  Fv  = -[bnu(  +  b12un]; 
G  =  J{F  •  Vt/)  =  G1  +  GV .  G1  =  a2u.  Gv  =  -[b2m(  +  b22u^}; 


~ai=Ja-  V£  bn  =  JV£  ■  BVi  bl2  =  JVi  ■  BVy 
d.2  =  Ja  *  Vtj  621  =  JVt]  -  BV£  622  =  JVr]  •  BVr\ 


and  B  is  given  in  (4.8). 

4.2.1.  The  energy  method.  Let 


(4.23) 

(u’v)J  =  f0  /i(ut;) 

II 

CN^ 

(4.24) 

(n=v)  =  fQ  J  (uv) 

d^drj. 

II 

CN 

IT 

(4.25) 

(n.  v)e.w  =  /  ( uv)e.w 

J  0 

dr). 

IMlI.w  =  0 u.u)E:w 

(4.26) 

(^t^)iV:5  =  J  {uv)N;S 

di. 

IMllr,s  =  {u-.u)N:s 

denote  the  weighted  L2  scalar  product  and  norm,  the  L2  scalar  product  and  norm,  the  boundary  scalar 
products,  and  boundary  norms  respectively.  The  subscripts  E.W.N.S  refer  to  the  EAST.WEST.NORTH 
and  SOUTH  boundaries  as  in  figure  4.1.  with  x.  y  replaced  by  £.77. 

The  equation  corresponding  to  (4.7)  becomes. 

(INIj)t  =  - .[(«,  F_ [  +  2 Fv)e  -  («,  FJ  +  2F^)w] 

EAST-WEST 

-  [(«.  G1  +  2 Gv)n  -  (ti,  G1  +  2 Gv)s] 

' - s, - ' 

NORTH-SOUTH 
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-  [K  -F/)  -  +  (U-  ^v)  ~  (u^:GJ)]+[(w^)+  (Mi)] 


GROWTH  1  GROWTH  2 

(4.27)  -[(UC;fy)+(>v!iiC)  +  KiGv)  +  (Gv;t/4. 

DISSIPATION 

Precisely  as  in  the  Cartesian  case.  GR1  and  GR2  in  (4.27)  can  lead  to  a  growth  or  decay  in  ||u||j.  but  will 
not  affect  well-posedness.  The  metric  relations  (4.19)  show  that  in  the  curvilinear  case  too..  GR1  vanishes 
for  constant  coefficient  problems.  Just  as  in  the  Cartesian  case,  we  need  to  assure  that  DI  has  the  right  sign. 

4.2.2.  The  dissipation.  Similar  conditions  as  in  the  Cartesian  case  for  positive  eigenvalues  also  apply 
in  the  curvilinear  case;  see  (4.27).  Thus  we  must  show  that 


bn  >  0;  622  >  0.  &11&22  — 


bi2  +  &21 


to  assure  the  right  sign  on  DI.  The  conditions  (4.28)  hold,  since  (4.22).  (4.8) .  and  (4.9)  lead  to  2 bn  = 
JV^t{B  +  Bt)V£  >  0. 2622  =  +  Bt)Vt)  >  0.  and 


611&22  ” 


&12  +  &21 


=  &11&22 


612  4-  &2i 
2 


4.2.3.  Boundary  conditions.  Consider  the  first  two  terms  in  (4.27)  and  recall  the  definitions  (4.3). 
The  outward  pointing  unit  normal  on  is. 

(4-29)  ntf  =  ±1,  „)  =  r,  =  0, 1)  =  ^ 

where  |V£|  =  yjil  +  Q  and  |Vjj|  =  +  Vy-  The  boundary  conditions  leading  to  an  energy  estimate 

become 

—ai(-l.Tj.t)  =  —Ja  ■  V£  <  0  =  F  —  Fw{r,,t).. 

—ai[—l.r].t)  =  —Ja-  V£  >  0  JFV  •  V£  =  Fv  =  F^fat), 

l'  '  +ai(+l.iM)  =  Ja  •  V£  >  0  JFV  •  V£  =  Fv  =  F^fat), 


-\-ai(-\-l.r}.t)  =  Ja-V£<  0  JF*V£  =  F  = 


at  £  =  ±1.  while 


-a2(£;0;t)  =  -Jff- Vr?<0  JF  -  V??  =  G  =  GsK,t): 

14  311  -fijK.0,1)  =  — Ja  •  V77  >  0  JFv-Vr,  =  Gv  = 

(  '  '  +a2(^.  l.t)  =  Ja  •  V»7  >  0  JFV -Vr)  =■  Gv  =  <%(£,*), 

+fi2(£>M)  =  Ja  •  <  0  JF-Vr)  =  G  =  Gat(£J) 

should  be  used  at  77  =  0. 1.  The  boundary  conditions  (4.30).  (4.31)  can  also  be  formulated  as  in  (4.12). 
The  same  procedure  as  in  the  Cartesian  case  leads  to  the  estimate 

(4.32)  (Hl5)r<  £  T||F,||?  +  GRl  +  GR2  +  DIi 

I=E.W.N.S  W 


Fe  =  o-i FE  +  (1  -  ai)Fl,  <n  =  -(1  -  |fii|/fii), 
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(4.33) 

and 


~  1 

Fw  —  osFw  +  (1  —  (T3)Fwi  0-3  =  —-(1  |ai|/ai). 

Gn  =  CT5Gjv  +  (1  —  05  )Gn:  0*5  =  -(1  —  |a2|/fl2): 
Gs  =  CTjG s  +  (1  —  07 )Gs  .  (J7  =  — 2^  1^2 1/«2) ; 


Ve 


fo 

Jo  u^dr) 


k=l: 


nw 


J01l5ll«2<ft? 

fo  u2drl 


/ix  |a2|«2«*€ 

VN  =  — Ti - T—  U=l,  ns-— i - - 1»=0- 

/_1«2<*£  /_  1«2^ 

The  parameters  7]e:  r}W:  TJN:  Vs  are  strictly  positive  if  a 1.  62  are  zero  for  a  finite  number  of  points. 

Time-integration  of  the  estimate  (4.32)  leads  to  an  energy  estimate  of  the  form  (2.3)  if  (4.28)  holds. 
Provided  that  a  solution  exists  we  can  conclude  that  the  following  theorem  holds. 

Theorem  4.  Problem  (4-20).  (4.30).  (4  31)  is  strongly  well  posed. 

4.3.  Interface  conditions.  Boundary  and  interface  conditions  of  the  flux  type  (see  (4.10) ;  (4.11). 
(4.30).  and  (4.31))  require  extra  careful  treatment;  see  [27]  for  an  example. 

4.3.1.  Interface  conditions  in  the  curvilinear  case.  To  apply  the  SAT  technique  [16]  on  the  fluxes 
at  an  interface  between  two  blocks  with  different  coordinate  transformations  and  matching  gridlines  (see  [17] . 
[18]  for  the  one-dimensional  treatment)  requires  that  we  identify  the  continuous  part.  Matching  gridlines  at 
£  —  £0  —  const  implies 

(4.34)  (®f)i  £  (x()2,  (y?)i  #  {y() 2,  (*„)i  =  (x„)2.  (yv)i  =  (ynh 

while  we  have 

(4.35)  (a;?)i  =  (xc)2:  (y$)  1  =  (y«)2:  (®„)i  ^  [xv)2.  (y„)i  ±  (y,)2 

at  t]  —  rjo  =  const.  The  subscripts  1.2  refer  to  the  two  coordinate  transformations. 

Equations  (4.21).  (4.19).  and  (4.34).  (4.35)  immediately  lead  to  the  conclusion  that 

(4.36)  Fid 0.t).t)  =  F2(£o.t).t).  Gida..r}.r)  ±  G2{Zo.t).t), 


(4.37)  Fid-.Vo-.T)  £  F2(tno:T):  Gid..T)0.T)  =  G2{t,rio,T): 

i.e..  F  is  continuous  across  £  =  const  while  G  is  continuous  across  77  =  const. 

4.3.2.  Interface  conditions  and  vanishing  wave  speeds.  Another  problem  with  flux-interface  con¬ 
ditions  appears  when  the  wave  speed  a  goes  to  zero.  Consider  the  two  constant  coefficient  problems 

ut  +  F(u)a;  =  0;  -L  <  X  <  0  and  vt  +  F(v)x  =  0.  0  <  x  <  L. 

where  F(w)  =  aw- h Fv(w).  Fv(w)  =  —  ewx.  Both  problems  have  homogeneous  outer  boundary  conditions 
at  \x\  =  L  and  zero  initial  data,  and  they  are  connected  through  interface  conditions  at  x  =  0.  We 
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will  compare  the  effects  of  flux-interface  conditions  (F(u)  =  F(v).Fv{u)  =  Fv (v))  and  variable-interface 
conditions  ( u  =  v.ux  =  on  the  solutions. 

By  transforming  the  problem  for  v  on  [O.+I]  onto  [-1.0]  via  the  transformation  x  -£.  and  then 
replacing  £  with  x:  we  obtain 


= 

Clpxx : 

t  >  0; 

-L  <  x  <  0. 

*  = 

0, 

t  =  0; 

—L  <  x  <  0. 

il 

1 

cq 

0, 

t  >  0. 

x  =  —L. 

B0tp  = 

0, 

t  >  0; 

x  =  0. 

where  ^  =  (w: r  A  =  diag(a;-a).  and  B-l$  =  0  denotes  the  outer  boundary  conditions.  Bq $  =  0 
represents  the  transformed  interface  conditions 

(4.39)  au  —  eux  =  av  -f  evxt  —eux  =  +evx  or  u  —  v.  ux~—vx. 

We  will  treat  (4.38)  as  a  half-plane  problem,  which  means  that  we  let  L  -*  oc  and  replace  the  influence  of 
B-l  by  only  admitting  bounded  solutions  as  x  — oc. 

The  Laplace-transfoxm  technique  applied  to  (4.38)  leads  to 


u(x.  s )  =  (T \ (s)  exp  (ki(s)z). 


v{x .  s )  =  02(5)  exP  (^2(5)2?) 


where  s  is  the  dual  variable  with  respect  to  time  and 


Kl~+Tt  + 


Note  that  both  u.  and  v  decay  away  from  the  boundary  x  =  0. 

The  interface  conditions  (4.39)  lead  to  the  equation  E[s)cr  —  0  where  <r  =  (^i:ct2)T-  A  well-posed 
bounded  solution  is  obtained  only  if  det(E(s))  /  0  for  5ft(s)  >  0.  The  flux-interface  conditions  in  (4.39  lead 
to 


(440)  £«)=(“;;;; 

while  the  variable- interface  conditions  leads  to 

(4.41)  E(s)  =  (  ^  dtHE[s))  =  2^/(^)2  +  7- 

Obviously  the  flux-interface  conditions  can  lead  to  unbounded  growth  for  vanishing  wave  speeds,  because 
det{E)a^  =  0  independent  of  s.  The  variable-interface  conditions,  on  the  other  hand,  lead  to  a  well-posed 
problem  since  det(E)a^ 0  =  2 y/Js/e). 

A  similar  analysis  of  the  flux-boundary  condition  an  —  eux  —  0  for  the  single  domain  yields  det(E(s))  = 
a/2~\~  y/(a/ 2)2  +  se.  Consequently,  the  problem  with  unbounded  growth  for  vanishing  wave  speed  does  not 
exist  in  the  boundary  condition  case  because  det(E)a_¥Q  =  y/ (sc). 

5.  The  discrete  problem  on  a  curvilinear  mesh.  In  the  rest  of  this  paper  we  will  consider  the 
transformed  problem  (4.20) ;  (4.30) .  (4.31).  Note  that  problem  (4.1).  (4.10) . (4.11)  corresponds  to  the  special 
case  where  r  =  t.  £  =  x  and  rj  =  y.  For  notational  simplicity,  we  ignore  the  “hat’;  notation  for  the  fluxes  and 
transformed  coefficients  introduced  in  (4.20)-(4.22). 
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NORTH 


£  =  -l  4  =  +l 

SOUTH 


FlG.  5.1.  The  single- domain  case  in  transformed  space . 

Let  the  iV  x  N  matrix  and  the  M  x  M  matrix  be  the  ID  positive  definite  matrices  defined  in 
Section  3.1.  A  product  av  is  arranged  discretely  (where  av  9a  AV)  in  the  following  way  (see  Figure  5.1): 


where  A{  —  diag(aij).  Also,  the  N  x  N  matrices  Je:  Jw:  and  the  M  x  M  matrices  Jjv;  Js:  Ip  are 


respectively.  The  subscripts  E .  W.  N.  S  refers  to  the  EAST.  WEST.  NORTH,  and  SOUTH  boundaries  (see 
figure  5.1).  At  this  point  we  also  define  the  restriction  of  U  to  the  boundaries: 

(5.4)  Ue:w  —  (Je.w  &  ir? )U:  Un.s  =  m  &  Jn;s)U. 

5.1.  The  norms  in  the  transformed  problem.  The  norms  and  scalar-products  corresponding  to 
(4.23)-(4.26)  are 

(5.5)  ( U.V)j=UT(P^Pv)JV,  (U;U)j  =  \\U\\2j, 

(5.6)  (U:V)  =  UT{Pt®Pv)V;  (U;U)  =  \\U\\2; 

(5.7)  [U.V)e:w  -Ut (Je.w  ®  Pi,)V  =  UeWPvVe,w ■  \\U\\e:W  =  {U.U)e:w : 
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NORTH 


£  =  -i  £  =  +i 

SOUTH 

Fig.  5.2.  Schematic  showing  J  requirements  necessary  for  MJ  to  be  a  norm:  r  —  4  and  q  =  2. 

(5.8)  (t/:  V>,s  =  UT(/\  ®  JN:S)V  =  UlsP(VN:S:  ||U||^-.5  =  (U,U)n:s- 

Obviously,  the  relations  (5. 6)- (5. 8)  define  norms  since  Pt £  and  Pn  are  positive  definite  matrices.  What  about 
(Pt®Pv)J  in  (5.5)  ? 

The  metric  scalar  J  is  defined  in  (4.19).  In  matrix  formulation  we  have 

(5.9)  J  =  diag(Ji).i  =  1. ....  N  Jj  =  diag(Jij),j  =  1. ....  M. 

We  need  the  following  theorem. 

THEOREM  5.  Let  M  —  P%  (&  Pv.  If  the  first  and  last  r  components  in  Ji  are  constants  and  the  first 
(J i ; ....  Jq)  and  last  Jn)  <1  blocks  in  J  are  equal  then  MJ  is  a  norm .  Proof:  Lemma  2.  with 

B  =  Pf:.  C  —  Prj  -  and  D  —  J.  with  J  defined  in  (5.9).  leads  directly  to  MJ  =  JM  =  M1/2JM1/2  >0  □. 

The  requirements  in  theorem  5  are  illustrated  in  figure  5.2  (for  r  =  4  and  q  =  2)  where  rxq  values,  of  J  are 
equal  in  the  corners  and  r.  q  values  of  J  are  constant  normal  to  an  rj  —  const.  £  —  const  boundary  respectively. 
A  curvilinear  transformation  with  such  a  J  close  to  JQ  is  called  volume  preserving  and  guarantees  that  MJ 
is  a  norm. 

Remark .  The  conditions  in  theorem  5  (i.e..  that  J  must  be  constant  in  the  first  q.  r  points  normal  and 
adjacent  to  the  boundary  JQ)  can  be  thought  of  as  theoretically  ideal  conditions.  In  practise  one  approaches 
the  ideal  condition  with  increasing  resolution  on  a  smooth  mesh  close  to  the  boundary  because 

j{ij)  -  J(0,j)  =  +  Q( A?):  i  =  1;  ...,9: 

+  0{Ai |2).  j  =  1. r. 

where  it  is  assumed  that  J(0.  j).  J(i.  0)  are  the  values  of  J  at  the  boundaries.  This  process  is  illustrated  in 
figure  5.3.  where  the  minimum  eigenvalue  of  PD  4-  DP  as  a  function  of  increasing  resolution  is  shown.  The 
minimum  eigenvalue  goes  from  a  negative  value  for  large  Az  to  a  positive  one  for  small  Ax. 

5.2.  The  single-domain  problem.  The  discrete  formulation  of  (4.20).  (4.30).  (4.31)  with  the  SAT 
technique  [16]  for  incorporating  flux  boundary  conditions  is 

(5.10)  JUT  4-  D^F  +  DVG  =  h  +  BC,  U( 0)  =  /, 

where  the  continuous  derivatives  F%.  Gv  are  approximated  with 

(5.11)  D(F  =  (P;1Qe<&In)Ft  DVG  =  (1%  &  P~xQn)G 
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GRID  REFINEMENT  ON  (PD  +DP) 


1 00  200 

Gridpolnts 


Fig.  5.3.  Minimum  eigenvalue  of  PD  -f  DP  as  a  function  of  Ax. 


BC  =  [Pf'Js  ®  i)(F  -  Ps)  +  (Pf1*  ®  I^2)(Fv  -  F%) 

+  (Pf 1  Jw  «  J„Ea)(P  -  JW)  +  (Pf1^  ®  J,E4)(Pv  -  P&) 

+  (JCEB  «  P^Jn)(G  -  Gn)  +  (7{E#  ®  P~1Jn)(Gv  -  G&) 

(5.12)  +(/^7®P-Vs)(G-G5)  +  (/{E8(gp-1Js)(GK-G^). 

The  N  x  N  matrix  and  the  M  x  M  matrix  Qv  are  defined  in  Section  3.1.  Fluxes  with  subscripts  E .  W.  AT.  S 
are  boundary  data.  The  matrices  Ei  -  Eg  will  be  determined  below. 

5.2.1.  The  energy  method.  Multiplying  (5.10)  from  the  left  with  PT(Pf  0  Pq).  introducing  the 
notation  M  =  P^  g  P^.  and  adding  the  transpose  of  the  equation  leads  to 

(||tf||5)r  =  -[t/T(Q*  «  PV)F  +  Ft(QJ  ®  P„)CT]  +  [{7TMft  +  /.TMi7] 

v - * - X  ' - * - ' 

A  GR2 

(5.13)  +  -[UT(Pt  «  Q„)G  +  GT(Pf  «  <#)17]  +BT  +  ( BT)t 


BT  =  Ut(Je  ®  P„Ei)(F  -  Fb)  +  t^T(/js  ®  Pr,E2)(Pv  -  Fb) 
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+  UT(JW  «  PVZ3)(F  -  Fw)  +  UT(JW  ®  PVU)(FV  -  Fw) 

+  UT(Pz E5  <8  Jn){G  -  Gn)  +  Ut{P&6  «  ^)(GV  -  G&) 

(5.14)  +£/t(PcE7®  J5)(G-Gs)  +  (7T(Pei:8«  Js)(Gv-G£). 

In  (5.13)  we  have  assumed  that  the  metric  transformation  is  such  that  MJ  is  a  norm.  Note  that  the  flux 
terms  with  subscripts  are  given  data. 

The  notations  and  abbreviations 

(5.15)  H-  Qj  =  —  Je  —  Jw :  Qtj  4-  —  Br]  —  Jn  —  Js t 

will  be  used  to  expand  the  A  and  B  in  (5.13).  We  obtain 

A+B  =  -  [U T{Bi  «  P,)(Ff  +  2 Fv)  +  ( F 1  +  2 Fvf  {B(  ®  P„)U]/ 2 
^  ^  v 

E-W 

-  [UT(P(  ®  Bf))(G/  +  2GV)  +  (GJ  +  2G'T(P?  ®  B,)tl]/2 

V  -  .✓ 

N-S 

-  {[(17,  D(FT)  +  [D^F1,  Cl)]  -  [( F D$U)  +  (D(U,  F1)}}/ 2 

' - V - / 

GR1 

-  {[(17,  D.G1)  +  (DvGr .  U )]  -  [(G1..  DnU)  +  (D„U:  G1)]}/ 2 

^  _  -  . > 

GR1 

(5.16)  +  [(D(U,FV)  +  {FV,D(U)  +  ( DVU).GV )  +  (GV1DVU)]. 

S  ^  .  ^ 

DI 

Note  the  close  similarity  of  the  discrete  energy  estimate  (5.13) ;  (5.14).  (5.16)  with  the  corresponding 
continuous  one;  see  (4.27).  Just  as  in  the  continuous  case  GR1  and  GR2  will  at  most  create  an  exponential 
time  growth.  To  obtain  an  energy  estimate  we  must  determine  under  what  conditions  the  dissipation  (DI)  is 
negative  definite  and  which  values  we  must  assign  to  the  matrices  £i  —  Y.&  to  obtain  bounded  contributions 
from  the  boundary. 


5.2.2.  The  numerical  dissipation.  The  DI  is 

(5.17)  DI  =  {D^U)tMFv  +  (FV)T  MD^U  +  (DnU)T  MGV  +  ( GV)TMDVU . 

The  relationship  between  the  gradients  and  the  fluxes  in  the  continuous  case  are  given  in  (4.21).  In  matrix 
formulation  for  the  discrete  case  we  have 

(5.18) 

where  Bki{i:j)  =  6/s/(&: rjj).  By  introducing  (5.18)  into  (5.17)  we  get 


■  FV  • 

’  Bu 

B\2 

'  D(U  ' 

Gv 

Bn 

B22 

.  DvU  . 

(5.19) 


DI  =  — 


D^U 


DnU 


B\\M  -j-  M  B\ i  ~b  M B\2 

B12M  +  MB21  B22M  -f  M  B22 
At  this  point  we  need  to  define  M  —  ®  P-q  in  greater  detail;  we  have 


I 

£ 

<=: 

_ 1 

.  DnU  \ 

(5.20) 


Pt  = 


1 

C 

<XJ> 

1 _ 

1 - 

c 

•if 

1 _ 

1  u 

1  u 

,Pn  = 

0 

0  1 

1 

_ 1 

— 1 

_ i 
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We  will  need  the  following  lemma. 

Lemma  3.  If  the  blocks  H%.Hf  have  the  size  q  x  q.  the  blocks  have  the  size  r  x  r.  and  the 

matrices  Bui  in  (5.18)  are  constant  in  the  first  q .  r  points  normal  and  adjacent  to  the  boundary  59l}  then  the 
dissipation  DI  defined  in  (5.19)  is  negative  definite.  Proof:  Lemma  2  leads  to 

B\\M  +  M  Bn  32i  M  +  MB\ 2 
B12M M B21  B22M  +  MB22 
M1/2  0  IT  2  Bn  B12  +  B21  1  r  M1/2  0  1  Q 

0  M1/2  B12  +  B2 1  2B22  J  [  0  M1/2 

Remark.  The  conditions  in  lemma  3  (i.e..  that  the  matrices  Bki  in  (5.18)  are  constant  in  the  first  q.r 
points  normal  and  adjacent  to  the  boundary  5Q )  can  be  thought  of  as  theoretically  ideal  conditions.  In 
practise  one  approaches  the  ideal  condition  with  increasing  resolution,  smooth  coefficients  bij  and  a  smooth 
mesh;  see  the  Remark  on  J  in  Section  5.1. 

5.2.3.  Stability.  To  obtain  an  energy  estimate  (given  the  negative  contribution  of  the  DI  on  the  right- 
hand-side  of  (5.16))  we  must  assign  values  to  the  matrices  Ei  -  E8  in  order  to  obtain  a  bounded  boundary 
contribution. 

Let  us  start  by  estimating  the  terms  at  the  EAST  boundary.  We  have 

BTe  =  -  {UT(PV(I/  2  -  Ei)]F7  +  (F7)T[(//2  -  I%)P,]U} 

-  { UT(Pn(I  -  Ei  -  Z2)]Fv  +  ( Fvf((I  -  E?  -  E  l)Pn]U) 

(5.21)  -  [UTPnFE  +  (Fe)tP,U]: 

where  Fe  =  Ei  Fe  +  E2  Fj . 

Obviously,  the  terms  involving  the  viscous  fluxes  must  be  removed.  This  yields  E2  =  /  -  £  1*  By 
observing  that  F 1  =  Aj %U  where  A&  =  diag((ai)jsjj)  (see  (4.21)  for  a  definition  of  a{)  we  obtain. 

BTe  =  -  UT[PV(I/ 2  -  Ei)]Ab  +  AE(//2  -  E J)PV]U 

(5.22)  -[VTPnFE  +  {FE)TP*U]. 

Now  we  choose  £1  such  that  (7/2  -  Ei)Ab  =  |AE|/2.  This  choice  and  an  entirely  similar  procedure  at  the 
other  boundaries  yields 

(5.23)  (||17|&)t<  £  +  GR1  +  GR2  4-  DI. 

I=E:W:N;S  ”J 

where 

FE  =  TiFE  +  (Iy-Ti)F%;  Ex  =  (Iy  —  |As|A^1)/2. 

Fw  —  E3  Fw  +  (Iy  —  E  3)7^;  E3  =  —(Iy  +  |Aiv|Aw/1)/2. 

Gn  =  E5G*  +  (Ix  -  E 5)G£:  E5  =  (Ix  -  |AJV|A^1)/2; 

(5.24)  Gs  =  E7Gs  +  (/,  -  E7)G£;  E7  =  -(/*  +  |As|Aj1)/2, 

(5.25)  E2  —  Iy  —  El .  E4  —  — Iy  ^3;  ^6  —  Ix  ^5;  ^8  —  f v  ^7; 

and 

l(U..MU)i  +  (\Ai\U..U)i]  I=EmNS 

m  -  2  (U.U)i  :  ;  ’  ’ 
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The  similarity  of  the  discrete  energy  estimate  (5.23)  with  the  corresponding  continuous  one  (see  (4.16) 
and  (4.32))  implies  strict  stability.  Time-integration  of  (5.23)  leads  to  an  estimate  of  the  form  (2.4)  if  the 
DI  has  the  right  sign.  i.e..  if  lemma  3  holds.  We  can  conclude  that  the  following  theorem  holds. 

Theorem  6.  The  approximation  (5.10)  of  the  problem  (4-20),  (4-30).  (4.31)  is  both  strictly  and  strongly 
stable  if  lemma  3  holds  and  Ei  —  Y.&  are  given  by  (5.24)  and  (5.25). 

5.3.  The  ID  multiple  domain  problem  revisited.  Before  we  consider  the  2D  multiple-domain 
problem,  let  us  once  more  look  at  the  ID  multiple-domain  problem  considered  in  [17].  [18]. 

5.3.1.  Derivation  of  the  Q-formulationfor  interface  problems.  Consider  the  following  hyperbolic 
interface  problem 


(5.26)  ut  +  Ux  =  0,  -1  <  x  <  0,  and  vt  +  vx  =  0.  0  <  x  <  1 

augmented  with  suitable  initial  and  boundary  data  and  the  interface  condition  u  =  v  at  x  =  0.  The 
straightforward  approximation  of  (5.26)  is 

Ut  +  P£1QlU  =  P£\<tl(Un  -  V0)eN) 

(5.27)  Vt  +  P^QrV  =  PrH'rIVo  -  UN)e 0) 

where  U  =  (Uo- ■■■■.  Un)t .  ejv  =  (0: ...,  0. 1)T,  V  =  (Vo; ...;  Vm)t o  =  (1»  0....  0)T. 

The  boundary  terms  from  the  left  (L)  and  right  (R)  outer  boundaries  are  ignored.  The  formulation 
(5.27)  can  also  be  written  in  the  following  way: 


(5.28) 


PWt  +  (Q  +  T)W  =  0 


where  W  =  ( U .  V)T .  P  =  diag(P f.  Pr):  Q  =  diag(QL.  Qr):  and 


(5.29) 


E  = 


0 

0 


0 

0 


-VL  +CL 
+(Tr  ct  r 


We  can  now  split  up  Q  +  E  into  a  symmetric  and  a  skew-symmetric  part  as 

„  .  „  (g  +  E)-(Q  +  E)T  ,  (Q  +  E)  +  (g  +  E)T 

y  T  ^  0  *  r> 


QB 


The  2x2  blocks  of  Qsk  and  D  corresponding  to  the  nonzero  elements  in  S  are 


1 

0 

l 

b 

1 

1  -  2  <7L  &L  +  CR 

2 

_  ~(crL  -  ur) 

0  J 

2 

<?l  +  a-R  -1  -  2ctr 

In  the  sequel,  the  “tilde”  sign  will  indicate  the  4  x  4  block  that  couples  the  solutions  in  the  left  and  right 
domains.  Equation  (5.28)  now  becomes 


(5.30) 


PWt  +  (Qsk  +  D)W  =  0. 


In  [17]  it  was  shown  that  (5.27)  is  conservative  if  ctr  =  ctl  “  1-  By  introducing  this  condition  in  Qsk  and  D 
we  obtain  the  final  form  of  the  difference  operator 


(5.31) 


O 

1 - 

1  ' 

l  -l  " 

.  b  —  a 

L-1 

o 

-l  i 
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where  a  =  1/2  —  07,. 

The  formulation  (5.30) .  (5.31)  hereafter  referred  to  as  the  Q-formulation  is  a  rearranged  form  of  the  orig¬ 
inal  formulation  (5.27).  However,,  the  Q-formulation  simplifies  and  even  extends  the  possibility  to  formulate 
suitable  penalty  terms  for  second  order  derivatives. 

5.3.2.  The  Q-formulation  for  advect ion- diffusion  interface  problems.  Consider 

(5.32)  ut  +  F{u)x  —  0.  -1  <  x  <  0.  and  vt  +  F{v)x  =  0.  0  <  x  <  1 

where  F(w)  =  a(x.t)w  —  ewx  augmented  with  suitable  initial,  boundary,  and  interface  conditions.  An 
approximation  of  (5.32)  using  the  Q-formulation  is 

(5.33)  PWt  +  ( Q3k)(AW )  -  e(Qsk  +  D2)P~1(Qsk  +  D»)  =  JW 

where  W  —  {U.  V)T  and  P  =  diag(PL:  Pr).  The  matrix  A  has  the  values  of  a(x{.t)  on  the  diagonal.  The 
operator  Qsk  is  defined  in  the  previous  Section,  and 


(5.34) 


i  —  1.2.3 


as  in  (5.31).  The  dissipation  Di  is  formulated  as  acting  on  W.  which  is  a  more  general  formulation  that 
includes  penalty  on  the  flux  [cr\  =  cra(O.t))  as  well  as  penalty  on  the  variables. 

We  can  now  prove 

Theorem  7.  The  approximation  (5.33),  (5.34)  of  the  problem  (5.32)  with  the  choices 


(5.35) 


<Ti  <0,  Cr2  =  0.  cr3  =  0 


is  conservative  and  stable. 

Proof:  The  energy  method  applied  on  (5.33)  leads  to 

\\Wft  =  (VW.AW)  -  {V[AW).W)-2c{VW,VW)~  WTB{AW  -2eVW)  +IT 

V - - - ^  ^ - - - /  s - - - x 

GR1  DI  BT 

where  VW  =  P~lQskW  and  the  interface  terms  IT  are  defined  as 


(5.36) 


IT_  W~T  ’  2Di+  2eD2P~1D3  e(D2  -  Ds)  "  W 

”  [  VW  J  0  L  e{D2  -Ds)  0  J  L  VW  J  0 


The  growth  (GRl)r  the  dissipation  (DI)  and  the  ordinary  boundary  terms  (BT)  match  the  terms  in  the 
corresponding  continuous  estimate  perfectly.  The  choices  (5.35)  makes  the  term  IT  maximally  negative 
definite  and  leads  to  stability.  The  approximation  (5.33)  can  now  can  be  written 


(5.37) 


PWt  +  Qsk{AW  -  eP^QskW)  =  DXW , 


which  leads  to  conservation.  □ 

It  remains  to  identify  the  penalty  terms  and  corresponding  interface  conditions  and  find  out  whether 
(5.37)  is  sufficiently  accurate.  The  equation 

PWt  +  Q(AW  -  eP_1QW)  -  (£>1  +  (Q  -  Qsk)A)W  +  e{QskP~xQsk  -  QP~1Q)W 
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NORTH 


«  =  -l  4  =  +l 

SOUTH 


Fig.  5.4.  The  multiple  domain  mesh  in  transformed  space . 


is  a  formulation  of  (5.37)  in  the  usual  penalty  form.  Obviously,  the  penalty  terms  (denoted  by  PT)  indicate 
that  the  one-sided  first  and  second  derivatives  are  replaced  by  first  and  second  derivatives  involving  the 
adjacent  domain,  and  that  dissipation  is  introduced  via  D\. 

By  introducing 


Q*k-Q  =  Ar  A  =  — - 


-(2(7!  +  a)  0 

0  (2d!  —a) 


we  get  PT  =  AAW+e(A(P~xQW)  +  QskAW):  this  result  shows  that  the  corresponding  interface  conditions 
are  u  —  v  and  ux  —  vx.  Also,  because 


A4>  =  (0. 0)T.  A (P~lQ<t>)  =  {0[Axl~\.  Axr£l).OiAxl~\,  Ax^))T 


where  <j>  is  a  smooth  function  and  r  is  the  order  of  the  approximation  at  inner  points,  the  approximation 

(5.37)  is  accurate  enough. 

5.4.  The  2D  multiple- domain  problem.  In  this  Section,  an  interface  at  £  =  0  with  matching 
gridlines  (see  Figure  5.4)  is  considered.  Matching  gridlines  implies  that  the  number  of  points  in  the  r) 
direction  (M)  is  the  same  on  both  sides  of  £  =  0.  We  will  also  assume  that  =  P^  =  Pv.  which  implies 
that  Q £  =  —  Q-q .  Note  that,  in  general,,  the  difference  operators  Df? .  D ^  can  be  different  in  the  left  and 

right  domains  and  that  A£l  ^  A£h  and  Nl  /  Nr. 

A  multiple-domain  Q-formulation  of  the  problem  (4.20) ;  (4.30).  (4.31)  is 

(5.38)  JWt  +  D\kF  +  DqG  =  M^{D  0  T.PV)W  +  h  +  BCr  W{ 0)  =  / 

where  W  =  (U.  V)T .  The  solutions  in  the  left  (L)  and  right  (R)  domains  are  denoted  respectively  by  U  and 
V.  and 


(5.39)  =M~1(Qlk  «P„),  DV  =  M~\P^QV). 

In  (5.38);  BC  denotes  the  boundary  conditions  in  (5.10)  at  the  NORTH.  EAST.  SOUTH.  WEST  boundaries 
in  penalty  form,  h  is  the  forcing  function,.  /  the  initial  data,  and  F  —  F1  4-  Fv .G  =  GJ  -f  Gv  the  fluxes 
where 


(5.40) 


i.eoF1  =  AXW.  Fv  =  -(BnDfW  +  BuDiW). 
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(5.41)  G*  =  AaW,  Gv  =  -(BnD\kW  +  B22DnW). 

The  remaining  definitions  and  notations  used  in  (5.38)  are  Q'lk  =  +  A. 


(5.42) 

M  — 

1 - 1 

°  ^ 

S  ° 

i _ i 

,  J  = 

'  Jl  o' 

0  JR 

.  = 

‘  Qf  o' 
o  Qf  . 

■  0 

0  ■ 

■  0 

0  1 

(5.43) 

A  = 

A 

,  D  = 

D 

0 

0 

0 

o  J 

(5.44) 


'  pt 

0 

•  A  =  -i 

’  1  -1  " 

,  D  = 

1  -1  ’ 

0 

1 

2 

1  -1 

-1  1 

The  matrix  coefficient  E  will  be  determined  by  stability  requirements. 


5.4.1.  Conservation.  The  Q-formulation  automatically  leads  to  conservation: 

Theorem  8.  The  approximation  (5.38)  of  (4-20),  (4.30),  (4-31)  is  conservative. 

Proof:  Let  h  —  BC  =  0.  Multiplying  (5.38)  with  the  integration  operator  M  where  <j>  is  smooth,  and 
observing  that  Qf  =  -{Q\h)T  +  B$,  Qv  =  -Q?  4*  Bn  {B^.BV  are  defined  in  (5.15))  leads  directly  to 

fMJWt  -  {D\k<t>)TMF  -  (. Dn4fMG  +  <j>T{B e  ®  P„)F  +  «  B„)G  =  0. 

The  approximation  (5.38)  is  conservative;  i.e..  it  reverses  the  process  of  differentiation  (second  and  third 
terms  above)  and  leaves  information  only  at  the  boundaries  (fourth  and  fifth  terms).  □ 

5.4.2.  Stability.  In  this  Section  we  will  prove  the  following  theorem. 

Theorem  9.  The  approximation  (5.38)  of  the  problem  (4-20),  (4-30),  (4-31)  is  both  strictly  and  strongly 
stable  if  theorem  6  holds  and  EP^  4-  P^E  <  0. 

Proof:  The  energy  method  applied  to  (5.38)  yields 

(5.45)  \\W\\2j)  =  GR1  +  GR2  +  DI+  BTE  +  BTW  +  BTN  +  BTS  +  IT 

where  it  is  assumed  that  MJ  is  a  norm;  the  requirements  are  given  in  theorem  5.  The  boundary  terms 
BTe  4-  BTw  4-  BTn  +  BTs  are  exactly  the  same  as  in  the  single  domain  case  (see  (5.24)).  while  the 
operator  in  GR1.  GR2.  and  DI  is  replaced  by  D%k  defined  in  (5.39).  Strict  and  strong  stability  of  (5.38) 
follows  if 


(5.46) 


IT  =  WtD  «  (EP„  +  PVE)W  <  0. 


Because  D  >  0.  we  need  EP^  4-  P ijE  <0.  □ 

Remark.  Because  Pv  >  0.  E  <  0  with  the  first  and  last  r  elements  in  E  being  constants  would  satisfy 
condition  (5.46). 

6.  Numerical  experiments.  In  the  calculations  below,  we  have  used  the  fourth-  and  sixth-order 
schemes  reported  in  [17]  in  space  and  a  five-stage  fourth-order  RK  scheme  [32]  in  time. 
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Fig.  6.1.  Instability  due  to  vanishing  wave  speed  and  flux  interface  conditions. 


6.1.  Vanishing  wave  speed.  For  problems  with  a  realistic  geometry,  one  will  frequently  encounter  zero 
wave  speed  somewhere  in  the  field  due  to  the  variation  of  the  metric  coefficients,  the  variable  coefficients,  or 
(for  nonlinear  problems)  the  solution.  This  difficulty  (see  Section  4.3.2)  particularly  severe  in  one  dimension, 
is  exemplified  in  the  calculation  of  Burger’s  equation  shown  in  Figure  6.1. 

The  instability  that  develops  close  to  zero  wave  speed  when  using  a  penalty  on  the  fluxes  at  the  interfaces 
is  evident.  With  interface  conditions  applied  on  the  variable  instead  of  the  fluxes,  the  instability  disappears. 
Also,  if  one  scales  the  problem  such  that  U  Varies  between  1  and  3  instead  of  0  and  2  one  can  use  flux 
interface  conditions  without  any  sign  of  instabilities.  This  anomalous  behavior  associated  with  a  vanishing 
wave  speed  occurs  with  other  numerical  schemes,  and  is  typically  suppressed  by  adding  dissipation  (e.g.  the 
“Entropy  fix”  used  with  Roe  solvers) . 

6.2.  Error  growth  due  to  varying  coefficients.  Consider  the  following  ID  test  problem. 

ut  +  Fc  =  o.  [x.y]efl.  t>o 

u  =  /.  [x.y]  E  D.  t  —  0 
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SYSTEM  STABILITY:  EIGENVALUES 
X=  [1 +(H)/(N-1)]'"ph* 


FIG.  6.2.  The  error  in  the  growth  rate  for  different  transformations , 

The  problem  (6.1)  is  1-periodic  in  y  and  has 

(6.3)  i*i  (0  .y;t)  =  au2(0  :y.t).  «2(1  :y:i)  =  Pui(l.y.t). 

as  boundary  conditions  in  the  ^-direction. 

By  introducing  a  2D  curvilinear  mesh  we  obtain 

Jiir  +  (*)€  +  (6),  =  :o  K,ij]en,  r>o 

1  j  II  =  /:  r  =  0 

where  F  =  J^F.G  =  and  Q  =  [f .  17]  E  [0.1]  x  [0. 1].  The  problem  (6.4)  has  the  same  boundary 

conditions  as  (6.1). 

6.2.1.  The  energy  growth  in  ID.  The  energy  growth  for  the  ID  (y  =  0.  rjx  =  0)  version  of  (6.1)-(6.4) 
with 

(6.5)  a  =  1-fez.  6  =  —1  -f  ez.  a  =  1.  (3  —  ^/(l  +  e)/(l  —  e) 

leads  to  ||i/||2  =  — e||ii||2.  The  growth  rate  —e/2  corresponds  to  a  single  eigenvalue  on  the  real  axis  in 
the  continuous  spectrum.  Figure  6.2  shows  the  error  in  the  sixth-order  numerical  approximation  of  this 
eigenvalue  for  different  transformations  (z  =  z(£)).  Figure  6.3  shows  the  convergence  (in  an  L2  sense)  of  the 
seven  eigenvalues  with  most  accurately  converged  real  parts.  The  convergence  rate  in  both  Figures  6.2  and 
6.3  is  at  least  6. 
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SYSTEM  STABILITY:  EIGENVALUES 

Seven  Moat 'Accurate  Eigenvalues 


Fig.  6.3.  The  error  in  the  growth  rate  for  varying  wave  speeds . 

Even  though  the  resolved  eigenvalues  (and  eigenvectors)  converge  at  the  theoretical  rate  (see  Figures  6.2- 
6.3).  there  are  unresolved  eigenvalues  and  eigenvectors  that  can  generate  difficulties.  In  Figure  6.4.  the  least 
resolved  eigenvector  corresponds  to  an  eigenvalue  with  a  negative  real  part  (-4.6529E-03)  significantly  more 
to  the  right  of  the  analytical  value  (-7.5000E-03)  than  could  be  expected  by  the  order  of  the  approximation. 
These  unresolved  eigenvalues  and  eigenvectors  may  generate  extra  large  energy  growth,  as  shown  in  Figure 
6.5.  The  growth  varies  with  the  initial  condition.  Note  that  the  extra  energy  growth  for  a  uniform  mesh 
can  be  present  only  for  varying  coefficients  because  otherwise  GR1  =  (F.  DXU )  —  {DXF.  U)  =  0. 

6.2.2.  The  energy  growth  in  2D.  The  energy  growth  for  the  2D  continuous  problems  (6.1).  (6.4) 
is  identically  zero  with  e  =  0  in  (6.5);  i.e..  the  L 2  norm  of  the  solution  remains  constant  in  time.  In 
the  semi-discrete  case,  the  energy  growth  is  given  by  (5.45)  where  GR2  =  DI  =  0  and  the  introduction  of 
boundary  conditions  ( BTi )  and  interface  conditions  (IT)  leads  to  damping.  Possible  error  growth,  see  (5.16). 
is  provided  by. 

(6.6)  GR1  =  -[(U.'DzF)  -  {D&F)]  -  [( U:DVG )  -  {DnU,G)\ 

only.  For  a  linear  mapping  where  the  metric  coefficients  are  constants  (see  Figure  6.6)  we  obtain  D^F  — 
D^AiU  =  Ai Df:U,  DVG  =  -  A2Dr7^;  which  yields  GR1  =  0.  The  error  growth  is  shown  in  Figure 

6.7.  The  calculations  are  fourth-order  accurate  in  time.  Note  that  there  is  an  absolute  bound  on  the  error. 

In  a  nonlinear  mapping  (see  Figure  6.8)  the  truncation  errors  in  the  metric  calculation,  and  consequently 
also  in  the  calculation  of  the  fluxes,  leads  to  GRl^  CL  which  in  turn  can  generate  error  growth  (see  Figure 
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Physical,  First,  and  Re<max) 


FlG.  6.4.  Eigenvectors  related  to  two  most  and  least  resolved  eigenvalues . 

6.9).  These  calculations  also  are  fourth-order  accurate  in  time.  Note  the  enormous  time  scale  in  Figures  6.7r 
and  6.9. 

6.3.  Navier-Stokes  calculations.  We  consider  here  a  ID  viscous  shock  propagating  in  accordance 
with  a  Mach  number  of  2.0  and  a  Reynolds  number  150  over  a  2D  domain.  The  exact  solution  of  the 
Navier-Stokes  equation  for  this  case  can  be  found  in  [33].  At  the  artificial  boundaries,  including  the  circular 
region  in  the  middle,  we  impose  flux  boundary  conditions  by  using  the  penalty  formulation  on  the  fluxes 
with  exact  data  from  the  analytical  solution.  At  the  interfaces  we  impose  interface  conditions  by  using  the 
penalty  formulation  on  the  variables. 

In  Figure  6.10.  the  density  and  grid  for  the  propagating  shock  is  shown.  The  shock  travels  from  the  lower 
left  corner  to  the  upper  right  corner  and  has  almost  passed  out  of  the  computational  domain  that  consists  of 
12  blocks.  The  sixth  order  scheme  and  24  gridpoints  were  used  in  each  sub-domain.  The  local  density  errors 
are  shown  in  Figure  6.11.  The  grid  refinement  study  in  Table  6.3  indicate  between  fifth-  and  sixth-order 
accuracy  in  an  L 2  norm,  consistent  with  the  theory  in  [34].  [35].  since  we  have  fifth  order  accuracy  at  the 
boundaries  and  interfaces  (see  (3.4))  and  relatively  coarse  grids. 

7.  Summary  and  conclusions.  We  have  analyzed  boundary  and  interface  conditions  for  high  or¬ 
der  finite  difference  methods  applied  to  multidimensional  linear  problems  in  curvilinear  coordinates.  The 
investigation  focused  on  the  effect  of  variable  coefficients. 

A  problem  with  a  norm  as  a  function  of  the  Jacobian  was  analyzed.  Boundary  and  interface  conditions 
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FlG.  6.5.  Extra  growth  due  to  unresolved  features  and  initial  conditions . 


Table  6.1 

Twelve  subdomains,  sixth- order  explicit;  CEL  =  0.3. 


Wave  speed 

49/65 

65/97 

97/129 

129/193 

-0.25 

-4.610 

-4.640 

-4.722 

-4.722 

0.00 

-5.115 

-4.986 

-4.538 

-4.657 

0.25 

-5.155 

-5.253 

-5.179 

-4.952 

0.50 

-5.331 

-5.401 

-5.467 

-5.327 

0.75 

-5.523 

-5.514 

-5.590 

-5.565 

1.00 

-5.635 

-5.622 

-5.659 

-5.719 

average 

-5.228 

-5.236 

-5.193 

-5.196 

in  both  flux  and  variable  formulations  have  been  investigated.  Flux  boundary  conditions  lead  to  energy 
estimates  whereas  flux  interface  conditions  lead  to  difficulties  if  the  wave  speed  approaches  zero. 

A  new  and  simplified  so  called  Q-formulation  of  the  penalty  method  was  derived  at  interfaces.  The 
Q-formulation  simplifies  and  extends  the  formulation  and  implementation  of  derivative  conditions  in  both 
one  and  two  dimensions  at  interfaces. 

It  was  shown  that  varying  coefficients  can  cause  unbounded  error  growth  via  the  truncation  errors 
even  though  the  boundary  and  interface  conditions  are  implemented  in  a  stable  and  dissipative  way.  The 
error  growth  may  be  large  due  to  unresolved  features  in  the  solution.  Numerical  calculations  confirmed  the 
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FlG.  6.6.  A  4  block  mesh,  linear  mapping . 


theoretical  conclusions. 
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Fig.  6.11.  Local  error  levels  for  propagating  viscous  shock . 


31 


